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NUMERICAL CALCULATIONS OF TWO-DIMENSIONAL, UNSTEADY 


TRANSONIC FLOWS WITH CIRCULATION 
Richard M. Beam and R. F. Warming 
Ames Research Center 


SUMMARY. 


The feasibility of obtaining two-dimensional, unsteady transonic aero- 
dynamic data by numerically integrating the Euler equations is investigated. 

An explicit, third-order-accurate, noncentered, finite-difference scheme is 
used to compute unsteady flows 'about airfoils. Solutions for lifting and 
nonlifting airfoils are presented and compared with subsonic linear theory. 

The applicability and efficiency of the numerical indicial function method are 
outlined. Numerically computed subsonic and transonic oscillatory aerodynamic 
coefficients are presented and compared with those obtained from subsonic 
linear theory and transonic wind-tunnel data. Proposed areas for future 
investigation are indicated. 


INTRODUCTION 


The transonic flight regime is especially susceptible to dynamic and 
aeroelastic instabilities because of (1) the maxima of the lift-slope curves 
and (2) the large time lag between surface motions and aerodynamic forces that 
occur for bodies moving at near, sonic speeds. These facts, coupled with the 
recent development of aircraft with transonic cruise speeds, have renewed 
interest in the solution of unsteady transonic aerodynamic problems. 

Unsteady subsonic and supersonic aerodynamic coefficients generally have 
been satisfactorily determined from linear theory (refs. 1-5). The transonic 
aerodynamic coefficients for airfoils can be obtained from linear theory if 
(refs. 6-8) 


e << 1 , kz « 1 , M e << 1, kM e << 1 

CO OO 

and if either 

I M - l| » e 2 / 3 

I 00 1 

(1 

k » e 2 / 3 

where k is the reduced frequency of oscillation; e, the thickness ratio; and 
M m , the free-stream Mach number. This linearization is developed from the 
small perturbation of a steady-state flow that is uniform. If the airfoil is 



not thin, the steady-state flow field will not be uniform but will depend on 
the nonlinear aerodynamic effects. There is still the possibility of a lin- 
earization of the small perturbation unsteady motion about the nonuniform 
steady state. The resulting equations, however, are still formidable since 
the nonlinear equations must be solved for the steady state and the partial 
differential equations with spatially variable coefficients for the unsteady 
perturbation. (This latter linearization is discussed in section IV.) 

Analytical solutions for the nonlinear transonic steady flow equations 
are very limited and generally are restricted to approximate solutions to the 
small perturbation potential equation. The complexity of the unsteady non- 
linear equations (low-frequency small perturbation or all frequency large per- 
turbation) or the unsteady linear equations with variable coefficients (small 
perturbation about nonuniform flow field) further restricts the development of 
analytical solutions. The development of large .computers and the refinement 
of numerical methods for solving difference equations have led to the solution 
of many heretofore unsolved aerodynamic problems. The numerical solution of 
steady-state transonic inviscid flows in both two and three space dimensions 
is well established (refs. 9-12). Similar numerical schemes seem the most 
promising for obtaining solutions to complex unsteady transonic flow problems. 

Time accurate numerical differencing schemes are available (refs. 9 and 
13) (although their primary application has been to solve steady-state equa- 
tions 1 ) and are probably the most efficient (least computer time) method for 
solving the unsteady nonlinear equations. While relaxation methods generally 
offer more efficient solution to steady-state transonic problems, they require 
one more dimension (artificial time) for time accuracy and would probably be 
less efficient than the time accurate methods. 

The most efficient method for numerically solving the linear variable 
coefficient equations (small perturbation about nonuniform steady state) is not 
evident. For time histories of general body motion, the time accurate methods 
probably are more efficient (for the same reason as in the nonlinear case) . 

If the desired result is the oscillatory aerodynamic coefficients, both time 
accurate solutions and relaxation methods have their advantages. The relaxa- 
tion methods do not require more dimensions than the time accurate methods 
since real time can be eliminated (harmonic motion) while artificial time is 
retained. Each frequency requires a separate solution to the relaxation 
problem. On the other hand, one solution to the time accurate method with an 
indicial (step) input for the motion can be used to obtain coefficients for 
all frequencies by a Fourier transform of the indicial response (which requires 
very little computer time). The time accurate indicial function approach is 
used here (sections IV and V) . The steady-state solution to the nonlinear 
equations, which provides variable coefficients for the linear problem, is 
obtained most efficiently by relaxation. However, the steady-state solution was 


^he limit of the time accurate methods as time becomes very large is the 
steady-state solution (ref. 9). For supersonic flows, the steady-state solu- 
tion is obtained more efficiently by taking advantage of the spatially hyper- 
bolic character of the equations (ref. 13). 
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obtained here from the same time accurate method (time approaching infinity) 
used to obtain the perturbation solution. 2 

This study is restricted to the two-dimensional flow of an inviscid per- 
fect gas. The Eulerian gasdynamic equations in conservation form were chosen 
rather than the potential equation. The choice of the form of the equations 
of motion was somewhat arbitrary and a similar study with the potential equa- 
tion seems warranted. The potential equation solution requires less computer 
storage and is generally the form best suited to relaxation methods. On the 
other hand, vorticity in the flow field (resulting from the unsteady motion of 
lifting bodies) and the embedded shock motion are probably most easily computed 
with the conservation form of the Eulerian equations. The most efficient 
choice (potential or Eulerian) remains to be determined. 

More approximate methods of solution (e.g., the mixing of supersonic and 
subsonic linear theory (refs. 14 and IS)) should not be excluded from consid- 
eration for transonic problems. The computation time required for such solu- 
tions could be considerably less than the present approach; however, consid- 
erable finesse is required in choosing the proper mixing. The solutions 
obtained. by the finite-difference methods should be useful in developing more 
efficient approximate methods. 

The numerical, explicit, third-order-accurate, finite-difference scheme 
and the equations of motion are given in section I. Previously unpublished 
details of the difference scheme are presented in the Appendix. Solutions for 
subsonic nonlifting and lifting airfoils are presented in section II and sec- 
tion III, respectively. Sections II and III indicate the adequacy of the 
numerical method in comparison with exact, linear, time-dependent solutions. 
Section IV reviews the use of indicial functions and their relation to the 
oscillatory aerodynamic coefficients. Numerically determined (via indicial 
function) oscillatory coefficients are compared with subsonic linear theory 
coefficients. Section V presents the results of a numerical analysis of a 
biconvex circular arc airfoil oscillating about midchord at transonic speeds. 
Comparisons are made between oscillatory moment coefficients obtained from 
linear theory, numerical solution, and wind-tunnel experiments. 


I. NUMERICAL METHOD 


The gasdynamic equations for two-dimensional, unsteady, inviscid flow 
are expressed in conservative (or divergence law) form as 





( 2 ) 


where t, x, and y represent the time and Cartesian space coordinates and U, F, 
and G are the vectors: 

2 In this study, the variable coefficient partial differential equations 
are not explicitly programmed. Rather the perturbation is made from the con- 
verged (nonuniform steady state) solution to the nonlinear equations. 
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( 3 ) 
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The variables in equations (3) are the density p, the velocity components u 
and. v, pressure p, and total energy per unit volume e. The additional 
required equation. is 'the equation of state (for a perfect gas): 

P = (Y - 1) [e - | (u 2 + y 2 )] (4) 

where y , a constant, is equal to the ratio of specific heats (taken as 1.4 
here). The conservative form of the Eulerian equations was chosen for its 
"shock-capturing" capabilities (ref. 16) . 

The numerical method is an explicit, third-order, noncentered, finite- 
difference scheme introduced by Warming, Kutler, and Lomax (ref. 17). The 
method is uniformally third-order-accurate in both time and spatial increments. 
The basic scheme applied to equation (2) takes the form: 



where 


U jlk = k^’kAy) , etc. (5d) 
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and w x and Wy are free parameters that may be selected to minimize dispersion 
or dissipation but must be properly selected to obtain a stable numerical solu- 
tion (see the Appendix and ref. 17.) 


During the development of the computer program, the' shock resolution was 
not satisfactory when the basic scheme (eqs. (5)) was used throughout the flow. 
In particular, a nonphysical oscillation near the shock produced a distorted 
pressure distribution upstream of the shock. Similar problems were encountered 
by investigators who used relaxation techniques (refs. 10 and 11). They found 
that the use of upwind difference schemes in the supersonic region resolved 
the difficulties. Accordingly, the program was modified to use the skewed 
upwind differencing scheme (described in the Appendix) for the supersonic por- 
tion of the flow. The method retains third-order accuracy in time and space. 

If j represents the streamwise (j increases downstream) direction, the skewed 
differencing scheme is based on points j-3, j-2, 3-1, j , and j+1. The specific 
differencing formulas are unchanged for the first and second steps; however, 
the final step (corresponding to eq. (5c)) becomes 
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where u xs is a free parameter similar to and oj 

For the present study, the flow tangency boundary condition is applied on 
the coordinate line instead of on the airfoil surface. This approximation is 
adequate for many steady transonic airfoil calculations. The boundary condi- 
tions for the far field were set equal to free stream. The effects of the 
outer boundaries were diminished and the resolution of the near-f^eld was 
increased by coordinate stretching. 

All computations were made on an IBM 360-67, coupled with an interactive 
graphics system. The graphics display was invaluable in checking out the 
program and in selecting the free parameters w x , and w xs and the optimal 
time step At. 
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II. FLOWS WITHOUT CIRCULATION (NONLIFTING AIRFOILS) 


M = 0.8 

BICONVEX CIRCULAR ARC 
(NONLIFTING) 


LINEAR THEORY 

o e « 0.02 
0 e = 0.06 


NUMERICAL 




To test the adequacy -of the numer- 
ical solution for subcritical noncircu- 
latory flows, the transient solution 
for a nonlifting biconvex circular arc 
airfoil was computed. The transient 
pressure coefficient distributions are 
presented in figure 1. The initial 
conditions were uniform free stream 
{M m - 0.8). The physical analogy is 
an airfoil accelerated from rest 
instantaneously to the free-stream 
velocity (step change in streamwise 
velocity component) . Numerical results 
are shown for two thickness-to-chord 
ratios (e = 0.02 and 0.06) at several 
time increments, t = U^t/c, where t is 
the number of chords traveled. The 
linear theory solutions were obtained 
by Lomax using methods similar to those 
in reference 3 and are presented for 
comparison. 

Several observations are useful 
in interpreting the results shown in | 
figure 1. The discontinuities of the ' 
pressure gradient (cp slope) predicted 
by linear theory do not appear in the 
numerical solutions for two reasons, 
one numerical and one physical. The 
numerical method tends to smooth these 
discontinuities because of the limited 
resolution (note that the linear theory 
pulse length is approximately the order 
of the numerical grid spacing) and 
numerical dissipation and dispersion. 
The physics introduces nonlinear terms 
in the equation of motion which are 

included in the numerical solution but are neglected by the linear theory. 

These nonlinea®ities can become important in regions of large pressure gradient 
changes, thus limiting the applicability of the linear theory. 

Note the significant difference in the two numerical solutions for dif- 
ferent thickness ratios. Note especially the lag of the receding wave in 
figures 1(b), (c) , and (d) and the greater pressure changes for the thicker 
airfoil. These effects should be anticipated because of the higher local Mach 
numbers and lower upstream propagation rates for thicker airfoils. The crit- 
ical pressure coefficient Cp (for M = 0.8) is 0.4346; thus the local Mach num- 
bers approach unity for the thicker airfoil (c?p/e = 7.24 for e = 0.06). 



x/c 



Figure 1.— Chordwise pressure coeffi- 
cient, o , at various times, fol- 
lowing a? indicial velocity change, 

U , at t = 0. 
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III. FLOWS WITH CIRCULATION (LIFTING AIRFOILS) 


One of the test cases for circu- 
latory flow was a flat plate with 
lateral velocity (sinking or plunging) . 
The initial conditions correspond to 
a step change in the lateraJL (normal 
to free stream) velocity component. 

The computed differential pressure dis- 
tributions for several time steps are 
shown in figure 2. The amplitude of 
the lateral velocity component for the 
numerical solution was chosen to cor- 



respond to an equivalent angle of 
attack (u/J/oo) of 1°. For this small 
perturbation, the agreement with 
linear theory should be good except 
for (1) the discontinuities in the 
pressure gradient predicted by the 
linear theory (as discussed in the 
previous section) and (2) the pressure 
coefficients near the leading edge 
where the linear theory predicts a 
weak singularity with infinite 
velocities . 

One advantage of using the full 
Eulerian equations (as opposed to the 
potential equation) is the ability to 
compute circulatory flow. For exam- 
ple, the transient motion of the sink- 
ing airfoil produces a vortex distri- 
bution in the wake that is swept 
downstream by the flow. Figure 3 com- 
pares the numerically computed wake 
circulation and that predicted by 
linear theory (ref. 3). The numerical 
values were obtained by evaluating a 
line integral around the wake at each 
time step of the numerical integration 
of the equations of motion. 

The section lift coefficient is 



obtained by integrating the chordwise 


pressure distribution. If the pressure 

distributions are those resulting from Figure 2.— Chordwise lifting pressure, 
a step change in velocity of the air- Ap, at various times, t, following 
foil, then the resulting lift coeffi- indicial sinking velocity, v ■, at 
cient time history C£ k (t) is called t = 0. 


the indicial lift coefficient. The 
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M = 0.8 

FLAT PLATE - SINKING 



’ T~ CHORDS TRAVELED 

Figure 3.— Circulation, r, fol- 
lowing an indicial sinking 
velocity, V , at t =0. 



Figure 4.— Indicial lift coeffi- 
cient, c ^ , for an indicial 

sinking velocity, v q , at x = 0. 


numerically determined indicial lift coeffi- 
cient for the sinking flat plate is shown 
in figure 4. The integrals of the numer- 
ically computed pressure distributions 
(fig. 2) were obtained by use of the trape- 
zoidal rule, while the linear theory inte- 
grals were obtained analytically and include 
the linear theory, leading-edge singularity. 

A discussion of the numerical oscilla- 
tions near x = 0 (fig. 4) is important since 
they indicate the maximum resolvable fre- 
quency of the numerical method. The fre- 
quency is proportional to the inverse of 
the grid spacing and the inverse of the 
Mach number. In terms of reduced frequency, 
k, the maximum resolvable frequency is of 
order 1/Af^A, that is, 

k ■ = q(w~t) (7) 

max WooV 

where A is the grid spacing made dimension- 
less with respect to the chord and M m is 
the free-stream Mach number. This does not 
present a severe limitation since a grid 
spacing of 1/20 the chord at M = 1 produces 
k max = 0(20), which is sufficiently above 
the range of most dynamic and aeroelastic 
investigations. The oscillations shown in 
figure 4 correspond to a reduced frequency 
of approximately 15. 


IV. INDICIAL FUNCTIONS AND OSCILLATORY AERODYNAMIC COEFFICIENTS 


Although one can integrate the nonlinear equations of motion (eq. (2)) 
for given initial conditions and time-dependent boundary conditions as in 
section II, the results have very restricted use since superposition of solu- 
tions (the basis of most flutter analysis) is no longer applicable. If oscil- 
latory solutions to the nonlinear equations are desired, the computation must 
be done for each frequency and each amplitude of oscillation and the coeffi- 
cients are valid only for the motion (or modal, function) assumed for the cal- 
culation (i.e., superposition of modal functions to obtain other motions is 
not applicable) . 
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The most productive area for the numerical investigation of unsteady 
transonic phenomena appears to be the small -perturbation unsteady motion about 
the nonuniform flow conditions (as outlined in the Introduction) (linear per- 
turbation about nonlinear steady state) . This type of linearization produces 
aerodynamic coefficients applicable to a particular airfoil geometry (as 
opposed to linear thin airfoil theory that has the same unsteady coefficients 
for all thickness and camber distributions) and the superposition of modal 
functions allows the application of conventional flutter methods. 

In the classical linear analysis of unsteady motion of airfoils, two dif- 
ferent but compatible approaches are taken. In one method, it is assumed ini- 
tially that the solution depends harmonically on time (refs. 1, 2, 4, and 5). 

For flutter analysis, the solutions (lift and moment coefficients) are tabu- 
lated as a function of reduced frequency and Mach number and are used in the 
airfoil equations of motion with the assumption that the motion is harmonic 
in time at the flutter (stability) boundary (ref. 18). If the solution is 
known for each frequency, the solution for the general time history is obtained 
with the aid of Fourier transforms (ref. 19). Although this method can be used 
to obtain general time histories, it is seldom used in practice since the 
oscillatory solutions are generally of primary importance. 

The second but less often used method for the linear analysis of unsteady 
motions is the indicial function approach (refs. 3 and 20). For example, if 
an airfoil is given an instantaneous sinking velocity (i.e., a discontinuous 
step change in velocity), the resulting flow field is the indicial response. 

For a given step change in motion (sinking, pitching, etc.), the solutions 
(generally in terms of lift and moment coefficients) can be tabulated as a 
function of time and, with the aid of Duhamel's integral, they can be used to 
compute the solution for a general time history of airfoil motion. The oscil- 
latory aerodynamic coefficients can be obtained by Fourier transforms of the 
indicial coefficients (refs. 3 and 20). 

Each approach can be used to obtain oscillatory aerodynamic coefficients 
by numerical solution of the equations for small perturbation about the steady- 
state nonuniform flow (linear partial differential equations with variable 
coefficients) . As mentioned in the Introduction, relaxation schemes appear 
most promising for solving the equations for harmonic motion (time eliminated) 
and time accurate methods appear best suited for the indicial function approach. 

For this study, the indicial function approach and time accurate solution 
of the Eulerian equations were exploited. The oscillatory aerodynamic coeffi- 
cients were determined by Fourier transform of the indicial response functions, 
which requires very little additional computational time. The indicial func- 
tions and oscillatory aerodynamic coefficients for sinking and pitching motion 
are related by: 
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The notation in equations (8) is that of Garrick and Rubinow (ref. 5) for the 
oscillatory coefficients (L\, L 2 > ^3> ^"+> ^l» ^2> M3, and M 4 ) and of Lomax, 
Fuller, and Sluder (ref. 3) for the indicial coefficients (c^, °m a > c ^q> and 

c m ) . Subscripts a and q refer to sinking and pitching motion, respectively. 
H 

The function Ac(x) is defined by Ac(x) = c(°°) - c(x). 


Figure 5 shows the results of a test case for obtaining the oscillatory 
coefficients using the numerically determined indicial functions. The indi- 
cial functions were computed for a sinking flat plate (see section III) and a 
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flat plate pitching about the lead- 
ing edge. These were used to com- 
pute the indicial function for a 
plate rotating about the quarter- 
chord point. This could have been 
obtained directly by computing the 
one indicial function for rotation 
about the quarter- chord point; how- 
ever, the two indicial functions 
allow computation of any motion in 
the plane. Equations (8) were used 
to obtain A/ 3 and M 4 plotted in fig- 
ure 5. Figure 5 also shows the sub- 
sonic linear theory results of Timman 
van de Vooren, and Greidanus 
(ref. 2), who solved the boundary 
value problem for harmonic motion 
using expansions in terms of Mathieu 
and modified Mathieu functions. The 
numerical results compare favorably 
with the results of Timman et at. 
and adequately predict the low- 
frequency destabilizing aerodynamic 
damping (positive A/ 4 , fig. 5(a)). 


The minimum frequency that can 
be resolved (k m £ n ) for a given amount 
of indicial function computation, 
T max> is not well defined and thus 
requires further investigation. 
However, the minimum frequency can be 
no less than 0(Tr/r max ) 


k . > 

min 


' T max ' 


( 9 ) 




Figure 5.— Oscillatory aerodynamic 
moment coefficients, A/3 and A/4, for 
various reduced frequencies k. 


The approximate range of resolution of reduced frequency can be stated from 
relations (7) and (9) as 


0 (r9 < * :<0 fe) (10) 

It was pointed out in section III that the step change in motion can 
introduce oscillations in the response curve (i.e., fig. 4). The oscillations 
are not physical but result from the introduction of a finite mesh size by the 
numerical method. They do not generally present an analysis problem since the 
frequency kmax (eq. (7)) must be kept greater than the maximum frequency of 
physical interest and the integral terms of equations (8) will effectively 
filter out the effect of the oscillations. If the amplitude of the oscilla- 
tions should become sufficiently large, the nonlinear effects in the equations 
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of motion could produce analysis problems. However, the amplitude of the 
oscillations can be minimized by the use of a ramp displacement function 
(finite rise time to steady value) in lieu of a step displacement motion. The 
rise time of the ramp motion x a (i.e., t q is the time to reach a constant 
amplitude of motion) should be chosen equal to the period of the induced oscil- 
lations (2 v/k max ) . The response function for the ramp motion and the response 
function for the step motion are related (with the aid of Duhamel's integral) 
by 


c R (t) = ~ f c S dx 1 , o < T < T 
T o " o ‘ l ' 

R (t) = f e 5 (Tl) dx l , T q < T < 00 


( 11 ) 


where superscripts R and S denote ramp and step, respectively. The frequency 
analysis using the ramp motion is now the same as for a step motion except it 
is generally more convenient to use the mean value approximation to evaluate 
the integrals in equations (8) rather than to differentiate (_ x) to obtain 

c^(t) . 


V. ANALYSIS OF CIRCULAR ARC AIRFOIL OSCILLATING 
IN TRANSONIC FLOW 


Subsonic linear theory predicts that an airfoil oscillating at low reduced 
frequency about a point at or forward of the quarter-chord point will have 
negative or destabilizing aerodynamic damping (ref. 4) (e.g., section IV and 
fig. 5(a)). For rotation about points aft of the quarter-chord, subsonic 
linear theory predicts stabilizing aerodynamic damping for all reduced fre- 
quencies. However, experimental wind-tunnel tests by Bratt and Chinneck 
(ref. 21) of a 7-1/2-percent thick biconvex airfoil oscillating about midchord 
exhibit an aerodynamic instability for a narrow range of subsonic free-stream 
Mach number (positive Wq, fig. 6(a)). The aerodynamic stiffness moment, m§, 
also decreases (fig. 6(b)), contrary to the prediction of linear subsonic 
theory. Experimental investigators believed that these effects were related 
to shock waves that formed on the surface of the airfoil; therefore, this 
model was chosen for numerical investigation. 

The procedure for obtaining numerical data points involved three steps at 
each Mach number investigated. First, the steady- state nonlifting solution 
was computed, then the indicial functions were computed, and finally the 
oscillatory .moment coefficients were computed using equations (8). Results of 
the numerical investigation are shown in figure 7. 

The experimental data of Bratt and Chinneck were obtained at very low 
reduced frequencies ( k = 0.004) — below the general range of flutter investi- 
gation and below the minimum reduced frequency ( k = 0.040) of the present 
numerical calculations (eq. (10)). To present a quantitative correlation 
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e = 0.075 





Figure 6.— Experimental oscillatory 
aerodynamic moment coefficients 
for biconvex circular arc airfoil 
rotating abount midchord. 


Figure 7.~ Oscillatory aerodynamic 
damping ( m •) and stiffness (m 0 ) 
moment coefficients for biconvex 
circular arc airfoil oscillating 
about midchord. 


between the numerical and experimental data obtained for different reduced 
frequencies, the subsonic linear theory results for the corresponding reduced 
frequencies are also shown in figure 7. The numerical results predict the 
negative damping and the decrease in stiffness moment obtained in the 
experiments . 
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Figure 8. — Chordwise oscillatory pressure 
coefficient and phase angle for biconvex 
circular arc airfoil oscillating normal 
to free stream (plunging). 


Chordwise pressure distri- 
butions were not measured by 
Bratt and Chinneck in the exper- 
iments discussed previously. 
Experimental pressure distribu- 
tions were obtained by Lessing, 
Troutman, and Menees (ref. 22) 
for a 5-percent-thick biconvex 
airfoil oscillating in bending 
(each spanwise section sinking) . 
Although the experimental model 
had a low aspect ratio {JR - 3.0) 
and there were three-dimensional 
effects, the chordwise pressure 
distribution at the quarter span 
station (midway between the 
fixed end and the tip) provides 
a qualitative comparison for 
experimental and numerical 
results (fig. 8) . The results 
are presented for the same 
reduced frequency and free- 
stream Mach number. The greatest 
disparity between the three- 
dimensional experimental data 
and the two-dimensional numerical 
results is, as anticipated, in 
the phase angle £, . While the 
qualitative agreement of pres- 
sure coefficient, , seems 

Pa h 

satisfactory, quantitative cor- 
relation must await three- 
dimensional numerical results. . 


CONCLUSIONS AND RECOMMENDATIONS FOR FUTURE INVESTIGATION 


Pressure distributions obtained by numerically integrating the two- 
dimensional unsteady Eulerian equations were presented. The applicability 
and efficiency of the numerical indicial function method were outlined. Sub- 
sonic and transonic oscillatory aerodynamic coefficients were computed and 
compared with those obtained from subsonic linear theory and transonic wind- 
tunnel data. The results of this study substantiate the feasibility of obtain- 
ing unsteady transonic aerodynamic data by numerically integrating the gas- 
dynamic equations of motion. 

The following topics are recommended for future investigations. . 
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(1) Since aerodynamic strip theory is not generally applicable in the 
transonic flow range, the present study should be extended to include three- 
dimensional flows. 

(2) The numerical integration of the unsteady potential equation should 
be investigated in a study similar to the present one. The lower computer 
storage requirement should be advantageous. 

(3) The optimum numerical method for obtaining oscillatory aerodynamic 
coefficients has not been resolved. Therefore, the application of relaxation 
methods should be investigated. The methods may prove particularly efficient 
for low reduced frequencies and for cases where coefficients are required at 
only a few reduced frequencies. 

(4) The more approximate methods for unsteady transonic flow analysis 
(e.g., those that mix the results of linear subsonic and supersonic theory) 
should also be pursued since they could require significantly less computer 
time. The results of the direct integration methods of the type used here 
should be useful in developing more approximate methods. 


Ames Research Center 

National Aeronautics and Space Administration 

Moffett Field, Calif. 94035, November 27, 1973 
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APPENDIX 


UPWIND THIRD-ORDER SCHEME 

This appendix describes an "upwind" third- order difference scheme that 
was used in the supersonic region of the transonic calculations described in 
this report. To review briefly, the basic noncentered third-order scheme was 
developed in reference 17. For a one-dimensional hyperbolic system in the 
conservation law form 


a u 
n 


+ 


wyj) = 0 

ax 


(Al) 


the scheme is 


UV -UP - 

3 .7 3 Ax 


n 


(A2a) 


= i (u. n ♦ u 

1 2 V 


) 


At Q 


( 1 ) 


Ax 


i/r l = a - 1 A* hi ? t 2 > . ■“ «■.(/” 

3 3 4 \ 3 / te j 4 Ax j 24 j 


(A2b) 

(A2c) 


where the conventional difference operators are defined by 


and 






v n - if 1 . 

3 3 - 1 









and 


The predictor L/^. 1 ^ and the first corrector are evaluated at time 

3 3 

(n+2/3)At and together constitute MacCormack's second order algorithm 
(ref. 17). The algorithm (A2) is called a noncentered scheme because it uses 
noncentered difference quotients to approximate aF/3x in the first two steps. 
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The fourth-order spatial operator <S 4 of equation (A2c) with the multiplicative 
parameter a> is necessary to achieve numerical stability. 

The von Neumann stability analysis of the scheme (A2) as described in 
reference 17 was based on the application of the scheme to the following 
linearized version of the hyperbolic system (Al): 


3 U . W 
at dx ~ 


o 


(A3) 


where A is the Jacobian matrix 3 F/W (A is assumed constant). A necessary 
and sufficient condition for stability is 


4v^ - v 


4 < 


go 1 3 


v| = \a\At/Ax - 1 


for all eigenvalues a of A. The range of go for a stable algorithm is shown by 
the shaded region in figure 9. The dashed curve of the figure is the value of 
w that minimizes dispersive error (ref. 17). 

The scheme (A2) uses numerical 
data from the five spatial points 

X j' X j± 1' 311(1 X j± 2 3t time leVel n t0 
advance the solution one time step. A 
third-order upwind scheme can be 
devised that uses data at four points, 
say Xj, Xj ±l , and ^-_ 2 , to advance the 

solution. However, to avoid undue com- 
plication in programming logic, the 
condition was imposed that the formulas 

for and U ^ (as given by 

J J 

eqs. (A2a) and (A2b)) be the same 
throughout the computational region. 

This led to the construction of the 
following upwind algorithm that uses 

data from the five points x^ x. ±1 , x._ 2 , and x._ % to advance the solution: 



Figure 9.— Stable range of free param- 
eter go and v for basic third-order 
scheme. 


= same as right side of equation (A2a) 


(A4a) 


V ^ 
3 


l/l +1 

3 


same as right side of equation (A2b) 


= U 


n At / 
’ " 4 \ 


1 + A 






n 


(A4b) 


3 

Ax 


24 g 


(A4c) 
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where E is the shift operator EU™ = ^j+\' ^ ast difference operator on 

the right of equation (A4c) is a fourth- order operator that is analogous to 
the free parameter terra on the right in equation (A2c). 

The remainder of this appendix outlines a linear stability analysis of 
the upwind scheme (A4) . The procedure is essentially the same as for the 
basic third-order scheme (A2) (see ref. 17 for a more detailed discussion). 

To perform the von Neumann (Fourier) stability analysis, the algorithm (A4) .is 
applied to a linearized equation (A3). A further simplification accrues if 
one applies the algorithm to the linear scalar equation 


3 u 3 u 

3t + ° 3x 


0 


(A5) 


rather than equation (A3) . Whatever stability bound results by application of 
the algorithm to equation (A5) is valid for the linear system (A3) if c is 
replaced by a (the eigenvalues of A). When applied to the linear equation (A5) , 
the algorithm (A4) can be reduced to a single-step difference operator if one 
inserts equation (A4a) into (A4b) and the result into (A4c) ; thus, 

U n+1 = U. n - jv(4 + A + v4 V 2 )Vu.” + y V 2 (l + y V)V\i&U 
0 0^ 5 Q Z Z <7 


- i v 3 (! + J V) V6 2 u^. n 


(JL) 

24 


£V 4 w , n 
0 


(A6) 


Yl 

where v = c At /Ax. To determine the amplification factor g(k ), each term u. 

Yl ^ 

in equation (A6) is replaced by the k- th Fourier component v ( k ) exp (ikjAx) 

Yl 

where v ( k ) denotes the k- th Fourier coefficient. The coefficients at n + 1 
and n are related by 

v n+ \k) = g{k)v n {k) 

The stability of the linear difference equation is ensured if \g{k) \ S 1 for 
all values of k. The square of the modulus of the amplification factor for 
equation (A6) can be expressed as 

\g(k~)\ 2 = 1 - I* 2 2 5(s) 

where 

z = sin 2 (0/2) , 0 =kAx 

and S(_z) is the following quadratic in z: 
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5 ( 3 ) = 3 [to + v 2 (v 2 + 2)] - 2(1 - v) 2 [ 3to + v(2v 3 + 4v 2 + 5v - 6 )]a 

+ [—to 2 - 2vto(4v 2 - 3v + 2) + 3v 2 (l - 2v) (2v 3 + v 2 - 4v + 4)]a 2 
A necessary and sufficient condition for stability is that 

5(a) ^ 0 for 0 ^ a S 1 

The stability constraint on the free parameter to is found by evaluating 5(a) 
in the small and large wave number limit, i.e., k = 0 and k = ir/Atc: 

5(0) = 3[to + v 2 (v 2 + 2)] * 0 

5(1) = [-to + 4v(l - v 2 )][io + 4v 3 - 4v + 3] >0 
Hence, there follows 


-v 2 (v 2 + 2) s to s 4v(l - v 2 ) 


The two boundary curves for to are 
plotted in figure 10 as a function of 
v. Clearly, equation (A6) is a stable 
algorithm for 0 <, v < 1 with to = 0. 
This is not the case for the basic 
third-order scheme (A2) , where a non- 
zero value of to is essential for sta- 
bility (see fig. 9). 


If the free parameter term in the 
last step of equation (A4) is retained, 
then a practical criterion for choosing 
the parameter to can be deduced from the 
modified equation (ref. 17) . The modi- 
fied equation represents the actual 
differential equation solved when a 
numerical solution is computed using a 
finite-difference equation. For the 
upwind algorithm (A6) , the modified equation is 



Figure 10.— Stable range of the free 
parameter to and v = 1 upwind 
third- order scheme. 


3 u 
3 1 


+ o 


du _ y' 

dX ~ ^ 
p=4 


r -i 3 ^ u 
PtP) ^ 


(A7) 
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where the first two coefficients are 


y(4) = - mt + v2(v2 + 2)] 

rr ^ Ax 5 ,, . v(4v 3 + 4v 2 + 9v - 61T 

u(5) - im - V >L“ * - 1 5 L \ 

The even-order spatial derivatives on the right side of equation (A7) corre- 
spond to dissipative effects and the odd-order spatial derivatives correspond 
to dispersive effects. The dispersive error of algorithm (A6) can be mini- 
mized by choosing to so that q(5) is zero. In this case, 

w = -v(4v 3 + 4v 2 + 9v - 6)/5 

A curve corresponding to this formula is plotted in figure 10. Past experi- 
ence (refs. 13 and 17) indicates that oscillations in the neighborhood of a 
discontinuity (shock) can be reduced if w is chosen for minimum dispersion. 

One of the primary advantages of noncentered finite-difference schemes 
(refs. 13 and 17) is that a scheme constructed for one spatial dimension can 
be directly generalized to two or more spatial dimensions. Computation of the 
supersonic region of the two-dimensional transonic flow described in this 
paper used the basic third-order scheme (A2) in the ^-direction and the 
upwind third-order scheme (A4) in the x-direction. The algorithm is given in 
detail by equations (5a), (5b), and (6) of the text. 

For the applications in this report, the to corresponding to minimum dis- 
persion was selected for the basic scheme (fig. 9). The same criteria could 
not be used for the upwind scheme because of the nature of the stability 
boundaries (fig. 10) . The minimum dispersion value of to that corresponds to 
the maximum value of |v| is no longer a sufficient condition for stability 
and care must be taken to ensure stability for the minimum value of |v| as 
well. For example, the one-dimensional Eulerian equations have the eigen- 
values u, u + e, and u - c. For supersonic flow, let u = a + Ac (A > 0); the 
eigenvalues then become (2 + A )c, (1 + A )c, and Ac. If the upwind scheme is 
applied throughout the supersonic zone (i.e., u > c) , then A = 0, 

a. = 0, |v| . =0, and w must be set equal to zero for numerical stability 

min 5 1 1 min n 

(fig. 10). However, if the upwind scheme is applied only in regions where the 
local velocity exceeds the local speed of sound by a fixed amount (i.e., 

A ^0), then Ivl . ^0 and nonzero values of to may be selected to more 

closely approximate the value of w for minimum dispersion. 
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